package gameEngine.agents.distrebution.accessories;


/////////////////////////////////////////////////////////////////////////////
//Copyright (c) 2004 CenterSpace Software, LLC                            //
//                                                                      //
//This code is free software under the Artistic license.                  //
//                                                                      //
//CenterSpace Software                                                    //
//260 SW Madison Avenue                                                   //
//Suite #112                                                              //
//Corvallis, Oregon, 97333                                                //
//USA                                                                     //
//http://www.centerspace.net                                              //
/////////////////////////////////////////////////////////////////////////////


/// <summary>
/// Class NormalDist represents the normal (Gaussian) probability distribution
/// with a specifed mean and variance.
/// </summary>
public class NormalDist
{

 private static  double MACHINE_EPSILON = 1.0e-12;
 
 private  static double[] a = 
 {
   2.2352520354606839287e00,1.6102823106855587881e02,1.0676894854603709582e03,
   1.8154981253343561249e04,6.5682337918207449113e-2
 };

 private static double[] b = 
 {
   4.7202581904688241870e01,9.7609855173777669322e02,1.0260932208618978205e04,
   4.5507789335026729956e04
 };
 private static double[] c = 
 {
   3.9894151208813466764e-1,8.8831497943883759412e00,9.3506656132177855979e01,
   5.9727027639480026226e02,2.4945375852903726711e03,6.8481904505362823326e03,
   1.1602651437647350124e04,9.8427148383839780218e03,1.0765576773720192317e-8
 };

 private static double[] d = 
 {
   2.2266688044328115691e01,2.3538790178262499861e02,1.5193775994075548050e03,
   6.4855582982667607550e03,1.8615571640885098091e04,3.4900952721145977266e04,
   3.8912003286093271411e04,1.9685429676859990727e04
 };
 private static double half = 0.5e0;

 private static double[] p = 
 {
   2.1589853405795699e-1,1.274011611602473639e-1,2.2235277870649807e-2,
   1.421619193227893466e-3,2.9112874951168792e-5,2.307344176494017303e-2
 };

 private static double one = 1.0e0;

 private static double[] q = 
 {
   1.28426009614491121e00,4.68238212480865118e-1,6.59881378689285515e-2,
   3.78239633202758244e-3,7.29751555083966205e-5
 };

 private static double sixten = 1.60e0;
 private static double sqrpi = 3.9894228040143267794e-1;
 private static double thrsh = 0.66291e0;
 private static double root32 = 5.656854248e0;
 private static double zero = 0.0e0;

 private static double OneOverRoot2Pi = 1.0/Math.sqrt(2*Math.PI);

 //endregion Static Variables 


 //region Instance Variables ----------------------------------------------

 //the mean of the distribution
 private double mean;
 
 //the mean of the Variance
 private double Variance;
 
 private final static double EPSILON = 0.000001;

 private double sigma;
 // Once we have the mean and variance, the following are constants
 // that are needed to evaluate the pdf.
 private double oneOverSigma;
 private double oneOverSigmaSqr;
 // c_ is just oneOverSigma_*OneOverRoot2Pi
 private double c_;

 //endregion Instance Variables 


 //region Constructors ----------------------------------------------------

 /// <summary>
 /// Constructs a NormalDist instance with the given mean and variance.
 /// </summary>
 /// <param name="mean">The mean of the density.</param>
 /// <param name="var">The variance of the density. Must be positive!</param>
 /// <exception cref="Exception">Thrown if the variance is less than or equal to zero.</exception>
 /// <remarks>The variance of the distribution is the standard deviation squared.</remarks>
 public NormalDist( double mean, double var )
 {
	 this.setMean(mean);
	 try {
		this.setVariance(var);
	} catch (Exception e) {
		// TODO Auto-generated catch block
		e.printStackTrace();
	}
 }




 /// <summary>
 /// Gets and sets the variance of the density.
 /// </summary>
 /// <exception cref="Exception">Thrown if the variance is less 
 /// than or equal to zero.</exception>
 /// <remarks>The variance of the density is the standard deviation squared.</remarks>
 public double getMean() {
	return mean;
}

public void setMean(double mean) {
	this.mean = mean;
}

public double getVariance() {
	 return( sigma * sigma ); //sigma SQR
}

public void setVariance(double variance) throws Exception {
    if ( variance <= 0.0 )
    {
      String msg = "Expected variance > 0 in NormalDistribution. Found variance = " + variance+"";
      throw new Exception( msg );
    }
    sigma = Math.sqrt( variance );
    oneOverSigma = 1.0 / sigma;
    oneOverSigmaSqr = oneOverSigma * oneOverSigma;
    c_ = oneOverSigma*OneOverRoot2Pi;
}


 /// <summary>
 /// Returns the probability density function evaluated at a given value.
 /// </summary>
 /// <param name="x">A position on the x-axis.</param>
 /// <returns>The probability density function evaluated at <c>x</c>.</returns>
 public double PDF( double x )
 {
   double y = (x - mean);
   double xMinusMuSqr = y*y;
   
   // c_ is a constant equal to one over sigma times one over square root of 2 PI
   return c_*Math.exp( -0.5*xMinusMuSqr*oneOverSigmaSqr );
 }

 /// <summary>
 /// Returns the cumulative density function evaluated at a given value.
 /// </summary>
 /// <param name="x">A position on the x-axis.</param>
 /// <returns>The cumulative density function evaluated at <c>x</c>.</returns>
 /// <remarks>The value of the cumulative density function at a point <c>x</c> is the
 /// probability that the value of a random variable having this normal density is
 /// less than or equal to <c>x</c>.
 /// </remarks>
 public double CDF( double x )
 {
   // This algorithm is ported from dcdflib:
   // Cody, W.D. (1993). "ALGORITHM 715: SPECFUN - A Portabel FORTRAN
   // Package of Special Function Routines and Test Drivers"
   // acm Transactions on Mathematical Software. 19, 22-32.
   int i;
   double del,temp,z,xden,xnum,y,xsq,min;
   double result, ccum;
   double arg = (x - mean) / sigma;

   min = EPSILON;
   z = arg;
   y = Math.abs(z);
   if(y <= thrsh) 
   {
     //
     // Evaluate  anorm  for  |X| <= 0.66291
     //
     xsq = zero;
     if(y > MACHINE_EPSILON) xsq = z*z;
     xnum = a[4]*xsq;
     xden = xsq;
     for(i=0; i<3; i++) 
     {
       xnum = (xnum+a[i])*xsq;
       xden = (xden+b[i])*xsq;
     }
     result = z*(xnum+a[3])/(xden+b[3]);
     temp = result;
     result = half+temp;
     ccum = half-temp;
   }

   //
   // Evaluate  anorm  for 0.66291 <= |X| <= sqrt(32)
   //
   else if(y <= root32) 
   {
     xnum = c[8]*y;
     xden = y;
     for(i=0; i<7; i++) 
     {
       xnum = (xnum+c[i])*y;
       xden = (xden+d[i])*y;
     }
     result = (xnum+c[7])/(xden+d[7]);
     xsq = Math.floor(y*sixten)/sixten;
     del = (y-xsq)*(y+xsq);
     result = Math.exp(-(xsq*xsq*half))*Math.exp(-(del*half))*result;
     ccum = one-result;
     if(z > zero) 
     {
       temp = result;
       result = ccum;
       ccum = temp;
     }
   }

   //
   // Evaluate  anorm  for |X| > sqrt(32)
   //
   else  
   {
     result = zero;
     xsq = one/(z*z);
     xnum = p[5]*xsq;
     xden = xsq;
     for(i=0; i<4; i++) 
     {
       xnum = (xnum+p[i])*xsq;
       xden = (xden+q[i])*xsq;
     }
     result = xsq*(xnum+p[4])/(xden+q[4]);
     result = (sqrpi-result)/y;
     xsq = Math.floor(z*sixten)/sixten;
     del = (z-xsq)*(z+xsq);
     result = Math.exp(-(xsq*xsq*half))*Math.exp(-(del*half))*result;
     ccum = one-result;
     if(z > zero) 
     {
       temp = result;
       result = ccum;
       ccum = temp;
     }
   }

   if(result < min) result = 0.0e0;
   //
   // Fix up for negative argument, erf, etc.
   //
   if(ccum < min) ccum = 0.0e0;

   return result;
 }

}
